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1.  Introduction 


The  nonlinear  behavior  of  turbulent,  high  Reynolds  number  (Re)  fluid  flow  con¬ 
tinues  to  defy  detailed,  theoretical  description  for  all  but  the  simplest  of  physical 
configurations.  Laboratory  experiments  in  wind  tunnels  or  water  channels,  while 
providing  a  controllable  environment,  have  limits  on  the  scales  of  flow  that  can  be 
reproduced.  Field  observations  rapidly  become  cost  prohibitive  with  increasing  in¬ 
strumentation  density,  cannot  be  simplified  by  removing  complicating  physics  (e.g., 
radiative  transfer),  and  lack  control  over  boundary  conditions.  With  the  increasing 
availability  of  supercomputing  resources,  numerical  simulations  of  turbulent  flow 
have  a  significant  role,  not  only  in  probing  the  behavior  of  turbulent  fluid  dynamics, 
but  also  in  providing  high-resolution  weather  forecasting  at  Army-relevant  scales 
in  the  atmospheric  boundary  layer  (ABL).  The  work  described  here  is  to  further  de¬ 
velop  vortex  filament  methods  to  investigate  and  simulate  the  unique  complexities 
of  the  ABL. 

Most  turbulent  fluid  models  in  use  today  are  associated  with  grid-based  numerical 
methods  that  solve  modified  forms  of  the  Navier-Stokes  equations.  Direct  numeri¬ 
cal  simulations  (DNS)  solve  the  Navier-Stokes  equations,  including  all  motions  as 
small  as  the  dissipation  subrange.  The  computational  costs  associated  with  cover¬ 
ing  a  large  domain  with  high  spatial  resolution  are  significant,  and  generally  require 
reducing  the  scale  separation  between  the  largest  energy  containing  range  and  the 
dissipation  range  (i.e.,  reducing  Re).  In  addition,  important  atmospheric  physics, 
such  as  radiative  transfer,  may  be  excluded  to  maintain  computational  tractability. 

Large  Eddy  Simulation  (LES),  as  a  general  strategy,  seeks  to  explicitly  simulate 
fluid  motion  for  a  range  of  scales  from  the  largest,  energy-containing  scales  just 
into  the  statistically  isotropic,  inertial  subrange.  The  effects  of  smaller-scale  mo¬ 
tions  are  modeled  via  a  subgrid-scale  (SGS  or  subfilter-scale,  [SFS])  parametriza- 
tion.  In  gridded  implementations  of  LES,  the  subgrid  scale  motions  are  modeled 
as  an  added  diffusion  several  orders  of  magnitude  larger  than  the  molecular  diffu¬ 
sion,  a  practice  that  is  equivalent  to  decreasing  the  effective  Re.  In  a  similar  vein, 
the  solution  of  numerical  schemes  on  computational  grids  adds  additional  unphys¬ 
ical  diffusion,  either  implicitly  by  the  scheme  itself  or  through  required  explicit 
smoothing  to  combat  numerical  instability.  Increased,  unphysical  diffusion  can  im¬ 
pact  the  generation  and  evolution  of  small-scale  vortical  structures  within  the  flow 
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that  should  be  explicitly  resolved.  Furthermore,  these  SGS  parametrization  gener¬ 
ally  preclude  the  modeling  of  energy  back-scatter,  in  which  small-scale  structures 
organize  into  larger  structures.  The  inability  to  properly  simulate  these  structures 
can  be  especially  evident  in  atmospheric  settings  with  strong  stable  stratification, 
where  intermittent,  turbulent  patches  and  small  vortices  produced  by  even  minor 
terrain  elements  are  commonly  observed.  The  SGS  or  SFS  parametrization  also  re¬ 
quires  additional  treatment  near  surfaces  to  account  for  dissipation  and  momentum 
and  heat  fluxes  near  the  surface.  In  atmospheric  models,  the  treatment  is  usually 
based  on  Monin-Obukhov  (M-O)  similarity;  however,  sloping  terrain  and  thermal 
stratification  violate  some  of  the  basic  assumptions  of  this  approach.  M-0  simi¬ 
larity,  often  with  ad  hoc  corrections,  is  commonly  used  despite  its  shortcomings 
because  of  the  lack  of  alternatives.  Quantifying  surface  fluxes  of  sloping  terrain 
under  thermal  stratification  at  full  Re  is  an  open  question  in  atmospheric  science. 

Vortex  methods  are  a  class  of  Lagrangian  and  semi-Lagrangian  computational  fluids 
models.  Rather  than  based  on  a  fixed  grid,  the  methods  simulate  fluid  flow  through 
mutually  interacting  and  convecting  vortex  particles.  The  vortex  blob  method,  which 
uses  a  radially  symmetric  vorticity  distribution  for  each  particle,  has  been  designed 
with  mainly  a  view  toward  simulating  flows  in  which  viscous  effects  are  impor¬ 
tant  everywhere.  To  accomplish  this,  the  mean  spacing  between  the  vortex  particles 
must  be  less  than  their  core  radius  in  order  to  converge  to  solutions  of  the  viscous 
Navier-Stokes  equations.  The  blob  method,  as  a  DNS,  may  be  well  suited  for  lower- 
Re  situations  or  for  laminar  flows  at  high  Re.  In  addition,  Yokota1  accelerated  the 
performance  of  an  implementation  of  the  vortex  blob  model  using  NVIDIA  graph¬ 
ics  processing  units  (GPUs)  to  74%  on  over  4000  GPUs  on  the  cluster  in  Tsubame, 
Japan.  The  core  numerical  scheme,  using  the  Fast  Multipole  Method  (FMM),  was 
packaged  and  released  as  ExaFMM,  an  open-source  CUDA  framework  to  acceler¬ 
ate  FMM  calculations.  For  turbulent  engineering  flows  at  high  Reynolds  numbers, 
the  techniques  used  for  accommodating  vortex  stretching  in  blob  methods  tend  to 
be  inadequate,  leading  to  instability.  In  contrast,  the  vortex  filament  scheme  in¬ 
cludes  an  efficient  and  effective  means  for  calculating  vortex  stretching,  which  is 
an  essential  requirement  for  modeling  turbulent  flows. 

A  series  of  technical  reports  and  papers  are  planned  that  describe  the  adaptation  of 
the  vortex  filament  method  to  ABL  flows,  especially  related  to  complex  terrain  and 
urban  domains.  This  first  report  gives  details  of  the  basic  method  and  presents  re- 


Approved  for  public  release;  distribution  is  unlimited. 


2 


suits  from  simple,  thermally  driven  flows.  Follow-on  reports  will  cover  the  method¬ 
ologies  for  simulating  complex  geometries  and  urban  environments  with  the  Vortex 
Filament  Method  (VFM),  and  will  provide  the  details  of  the  FMM  as  it  has  been 
adapted  from  the  ExaFMM.2 

2.  Vortex  Filament  Method 
2.1  Overview 

The  VFM  solves  the  incompressible  vorticity  equation3 : 


(1) 


where  u>,  u,  p,  and  p  are  the  vorticity,  velocity,  pressure,  and  density,  respectively, 
The  left-hand  side  is  the  material  derivative  of  vorticity  and  the  right-hand  side  are 
the  vorticity  production  terms  that  generate  vorticity  through  stretching,  viscosity, 
and  baroclinic  effects.  Because  it  solves  the  vorticity  equations,  the  VFM  it  well 
suited  for  flows  dominated  by  vorticity. 

The  VFM  is  a  Lagrangian  method  where  the  computational  elements  are  small, 
straight-line  vortex  tubes  that  connect  end-to-end  to  form  vortex  filaments.  An  ex¬ 
ample  filament  is  shown  in  Fig.  1 .  Tube  end  points  provide  a  natural  discretization 
of  the  filament.  Each  filament  has  a  constant  circulation  that  is  determined  during 
creation  through  initialization,  boundary  conditions,  or  one  of  the  vorticity  pro¬ 
duction  terms.  This  circulation  induces  a  velocity  field  that  convects  the  tube  end 
points,  causing  the  tubes  to  stretch  or  shrink  and  the  filaments  to  fold.  Any  tube  that 
stretches  beyond  a  maximum  length,  lmax,  is  split  into  smaller  tubes  that  satisfy  this 
constraint.  Similarly,  if  a  tube’s  length  becomes  less  than  a  minimum  length,  lmin, 
it  is  removed  and  the  filament  is  reconnected.  Thus,  production  of  vorticity  through 
stretching  is  accounted  for  intrinsically  as  the  tubes  move  in  response  to  the  flow. 
The  process  of  stretching  and  folding  is  the  mechanism  by  which  energy  travels 
from  large  scales  to  small  scales. 
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Fig.  1  A  filament  consisting  of  small,  straight-line  vortex  tubes  connected  end-to-end 


The  VFM  is  adapted  to  thermally  driven  flows  by  solving  the  energy  equation  using 
a  distribution  of  energy  particles  that  represent  the  system’s  energy.3  Energy  and 
vorticity  are  coupled  through  the  velocity  field  and  baroclinic  vorticity  generation, 
which,  when  using  the  Boussinesq  approximation,  is  proportional  to  the  horizontal 
temperature  gradient.  At  each  time  step,  temperature  gradients  are  computed  from 
the  energy  particle  distribution  and  the  associated  vorticity  is  released  into  the  flow 
field  as  filaments. 

In  general  the  velocity  field  is  composed  of  the  velocity  induced  by  the  filaments, 
and  depending  on  the  problem,  a  free  stream  velocity  and  a  potential  flow  velocity 
for  imposing  a  no  penetration  boundary  condition.  For  the  work  presented  here, 
only  velocity  induced  by  the  filaments  is  considered. 

2.2  Velocity  Field 

The  velocity  field  induced  by  a  collection  of  vortex  tubes  is  governed  by  the  Biot- 
Savart  law  and  is  computed  using  the  information  carried  by  each  tube,  namely 
circulation,  length,  position,  and  orientation.  This  procedure  is  described  after  re¬ 
viewing  the  concepts  of  vortex  lines,  tubes,  and  filaments. 

A  vortex  line  is  a  curve  that  is  everywhere  parallel  to  the  vorticity.  Given  some 
closed  curve  C,  the  collection  of  all  vortex  lines  passing  through  each  point  of  C 
defines  the  surface  of  a  vortex  tube.  Circulation  T  of  the  tube  is  defined  as  the 
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integral  of  velocity  around  the  closed  curve  C  and  is  related  to  the  vorticity  by 


T  =  j)  u  ■  ds  = 

C  S  8 

where  vorticity  uj  =  V  x  u,  ds  is  the  differential  length  along  curve  C,  S  is  any 
surface  whose  boundary  is  C,  and  n  is  the  normal  to  S.  Since  V  •  cu  =  0  and  u>  is 
everywhere  parallel  to  the  surface  of  the  tube,  the  divergence  theorem  implies  that 
the  circulation  T  around  any  vortex  tube  is  a  constant.  The  theorem  of  Helmholtz 
and  Kelvin  shows  that  this  is  true  for  all  time.  A  vortex  filament  is  defined  by  a 
limiting  process  where  the  cross-sectional  area  of  a  tube  goes  to  zero.  In  this  limit, 
T  remains  constant  as  dS  — >•  0  and  I  uj  I— »  oo. 


V  x  u  ■  n  dS  =  u:  ■  n  dS, 


(2) 


Continuous  filaments  have  discrete  counterparts  in  the  VFM.  In  the  VFM  a  fila¬ 
ment  consists  of  small,  straight-line  vortex  tubes  connected  end-to-end  to  form  the 
filament.  A  maximum  length  is  imposed  on  the  tubes  so  that  a  discrete  filament  is 
allowed  to  bend  and  fold  as  an  approximation  to  a  continuous  filament.  Each  dis¬ 
crete  filament  carries  with  it  a  constant  circulation,  and  the  velocity  induced  by  it  is 
computed  using  a  distribution  of  point  sources  located  at  the  centers  of  the  tubes. 


Using  the  Helmholtz  decomposition  theorem,  Bernard3  shows  how  knowledge  of 
dilatation  and  vorticity  can  be  used  to  reconstruct  any  smoothly  varying  velocity 
field  in  an  unbounded  domain.  For  this  work,  incompressible  flow  is  assumed  so 
knowledge  of  vorticity  is  sufficient  to  compute  the  velocity  field.  The  formula  for 
this  is 


u(r,t) 


u>(s,  t)  x  (r  —  s) 
I  r  —  s  I3 


dV(  s). 


(3) 


R3 


Integrating  over  the  support  of  the  vorticity  in  R3  yields  the  velocity  at  r  at  time  t. 


In  the  VFM,  the  support  of  vorticity  is  captured  by  the  filaments.  Tubes  in  each 
discrete  filament  represent  a  small  volume  of  constant  vorticity.  Assuming  r  ^  s, 
integral  Eq.  3  is  used  to  approximate  the  velocity  contribution  at  r  of  a  vortex  tube, 
with  infinitesimally  small  volume  Sts  and  centered  at  s,  by 


Su(  r,  t) 


1  u>(s ,t)  x  (r-s) 

A  I  |Q 

4tt  r  —  s  d 


(4) 


The  point  source  approximation  is  obtained  in  the  limit  as  the  tube  volume  shrinks 
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to  point  s.  Because  circulation  T  is  constant  and  remains  so  in  the  limit,  rewriting 
Eq.  4  in  terms  of  the  tube’s  circulation  yields  the  point  source  formula.  Assuming 
the  tube’s  vector  length  is  dl  and  cross-sectional  area  is  dS,  then  using  the  relations 
Sts  =  dS  |  dl  |,  u)  =|  u)  |  dl/  |  dl  |,  and  r  =|  u>  |  dS  =  constant,  the  induced 
velocity  at  r  of  the  source  at  s  is  then 


Su(r,  t ) 


rdlx(r-s) 

47 r  I  r  —  s  I3  ’ 


(5) 


again  assuming  r/s.  The  influence  of  the  entire  filament  is  the  sum  of  all  the  tubes 
in  the  filament. 


2.3  Velocity  Smoothing _ 

A  vortex  core  is  characterized  by  a  circular  region  of  finite  radius  where  the  vor- 
ticity  is  constant  and  the  core  rotates  as  a  rigid  body.  If  the  numerical  method  used 
filaments  that  have  a  small  but  finite  core  radius  and  used  the  3-D  volume  integrals 
to  compute  the  velocities,  then  the  numerical  difficulties  associated  with  the  singu¬ 
larity  in  the  point  source  formulas  would  vanish  and  the  core  region  would  rotate  as 
expected;  however,  because  of  the  computational  cost  of  the  3-D  volume  integrals, 
this  is  not  practical.  Consequently,  the  core  region  is  typically  modeled  by  applying 
a  smoothing  function  to  Eq.  5,  see  for  example  Bernard4  and  Beale.5  If  <j)a  (|  r  —  s  |) 
is  the  smoothing  function  with  parameter  a,  then  Eq.  5  becomes 

t)  =  j- dl ,x  (r  |3S)  <Mlr-s|).  (6) 

47T  |  r  —  s  p 

The  smoothing  function  <pa  must  be  at  least  0(|r  —  s|2)to  eliminate  the  singularity. 
It  must  equal  1  or  quickly  approach  1  in  some  neighborhood  away  from  s,  and  allow 
the  velocities  to  go  to  0  as  |  r  —  s  |— »  0. 


It  is  interesting  to  compute  the  velocities  induced  by  a  single  vortex  tube,  as  used  in 
the  VFM  but  with  finite  volume,  using  the  3-D  volume  integral  Eq.  3  and  compare 
them  with  to  the  velocities  computed  using  the  smoothed  point  source  formula 
Eq.  6.  For  this  comparison,  2  smoothing  functions  are  used:  one  from  Bernard4 
given  by 


^(lr-sl)  =  1 


(7) 
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and  a  simple  piecewise  smoothing  function  given  by 


0a( I  r  —  S  I) 


if  |  r  —  s  |  >  a 
if  I  r  —  s  |<  cr. 


(8) 


Each  is  O  (|  r  —  s  |3)  and  satisfies  the  requirements  of  the  smoothing  function. 


Consider  a  single  vortex  tube  as  shown  in  Fig.  2,  which  is  centered  at  the  origin, 
aligned  with  the  a>axis,  and  has  core  radius  rc  and  length  l.  Let  the  vorticity  asso¬ 
ciated  with  this  tube  be  c o  =  (u>x,  0,  0)  and  let  x0  =  (x0,  yo,  0)T  be  any  arbitrary 
evaluation  point  in  the  z  =  0  plane.  In  this  plane,  the  only  nonzero  induced  velocity 
is  itz(x0)  and  is  given  by 


tc(x0) 


rc  1 /2 

^ x 

Alt 

~rc  -1/2 


_ (yo  -  y) _ 

[Oo  ~  X )2  +  (yo  -  y)2  +  z2]% 


dy  dx  dz. 


(9) 


The  first  2  integrals  are  integrated  exactly  and  can  be  found  in  the  Appendix.  A 
nonsingular  integral  remains  that  is  numerically  integrated  for  the  comparison  to 
the  smoothed  point-formula  velocities. 


Fig.  2  A  vortex  tube  centered  at  the  origin  and  aligned  with  the  x-axis  with  core  radius 
rc  =  0.001  and  length  1  =  0.0025 


Figure  3  compares  uz  computed  with  the  3-D  volume  integral  to  uz  from  the  point 
source  formula  with  both  smoothing  functions  (f^  and  <\}a  using  the  vortex  tube 
given  in  Fig.  2  with  circulation  =  0.001.  The  ^-velocity  is  plotted  as  a  function  of 
r/a  (r  =|  r  —  s  |)  for  2  different  approach  angles:  Fig.  3a  approaches  the  center  of 
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the  tube  along  the  y-axis,  and  Fig.  3b  approaches  it  at  a  45°  angle  from  the  positive 
x—y  quadrant,  both  in  the  z  =  0  plane.  As  expected,  beyond  2  core  radii  all  3  curves 
come  together  and  the  velocity  approaches  zero  as  r  gets  small.  The  peak  velocities 
near  the  core  radius  are  considerably  different  for  all  3,  the  largest  being  about  a 
116%  difference  between  (pba  and  the  volume  integral  velocities.  Both  smoothing 
functions  give  the  desired  behavior  but  if  compared  with  the  3-D  integral  results, 
the  linear  smoothing  function  seems  a  bit  better.  One  could  argue,  however,  that 
each  models  the  desired  behavior  correctly  but  is  based  instead  on  some  effective, 
but  different,  core  radius.  It  is  difficult  to  say  if  one  is  actually  better  than  the  other 
since  the  core  radius  is  a  scheme  parameter,  and  in  a  sense  so  is  the  smoothing 
function.  The  differences  may  have  some  effect  on  the  solution,  possibly  on  the 
Reynolds  number,  and  uncovering  this  relationship  is  one  of  the  long-term  goals  of 
this  research. 


Fig.  3  Comparison  of  uz  computed  using  the  3-D  volume  integral  (black)  with  the  point  source 
formulas  using  smoothing  functions  (red)  and  (.%  (green).  The  tube  configuration  is  given 
in  Fig.  2.  Velocity  is  plotted  as  a  function  of  y0  for  various  values  of  x0.  The  circulation  used 
was  0.001. 


2.4  Loop  Removal 

Simulations  of  complex  turbulent  flows  produce  convoluted,  spatially  intermittent 
vortical  structures  of  intricate  detail  that  tend  to  require  a  phenomenally  large  num¬ 
ber  of  elements  for  their  description  through  time.  The  vortex  stretching  process, 
which  is  accompanied  by  folding  that  brings  energy  to  small  dissipative  scales, 
leads  to  an  exponential  growth  rate  in  the  number  of  tubes  required.  Controlling  this 
growth  rate  is  absolutely  essential  for  any  simulation.  Bernard4  has  demonstrated 
that  loop  removal  reduces  this  to  a  linear  growth  rate  without  apparent  harm  to  the 
underlying  physics,  the  justification  being  that  direct  elimination  of  folded  vortices 
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in  the  form  of  loops  removes  primarily  local  energy  that  is  likely  destined  for  subse¬ 
quent  dissipation  at  smaller  scales.  In  his  conclusions,  Bernard4  acknowledges  this 
loss  of  energy  due  to  loop  removal  needs  to  be  further  investigated  to  uncover  the 
nature  of  this  relationship.  Uncovering  this  relationship  is  another  long-term  goal 
of  this  research  effort,  although  not  part  of  this  first  report.  The  expectation  is  that 
the  code  developed  will  eventually  be  capable  of  simulations  with  1  to  2  orders  of 
magnitude  more  vortex  tubes  than  previous  calculations  in  hopes  of  uncovering  this 
relationship,  and  others  associated  with  turbulent  flows.  The  actual  loop  removal 
algorithm  used  in  this  study  is  described  in  Algorithm  3  of  Section  3. 

2.5  Energy _ 

Adapting  the  VFM  to  thermally  driven  flows,  such  as  found  in  the  atmosphere,  re¬ 
quires  solving  the  coupled  energy  and  momentum  equations.  The  energy  equation, 
under  the  assumptions  of  constant  thermal  conductivity  and  no  deformation  work, 
reduces  to  a  convection-diffusion  equation  for  temperature3 : 

QT 

—  +  u  ■  vr  =  «v2t,  (io) 

at 

where  a  is  the  thermal  diffusivity. 


The  Boussinesq  approximation,  which  is  appropriate  for  many  atmospheric  flows, 
simplifies  the  momentum  equation.  It  ignores  variations  in  density  except  where  it 
is  multiplied  by  gravity.  Assuming  an  ideal  gas  undergoing  an  isobaric  process,  a 
first-order  approximation  to  density  can  be  written  as 

P  ~  Po  +  AT  =  po  (1  -  AT) ,  (11) 

where  AT  =  T  —  Tref  and  aVo  is  the  isobaric  coefficient  of  thermal  expansion 


1  (dp\  _  1  _  1  /dV\ 

~p[df)p~f-v[df)p 


(12) 


evaluated  at  (po,po).  Ignoring  variations  in  density  in  all  but  the  gravity  term,  the 
momentum  equation 


P 


+  (Vu)  u 


pg  -  Vp 


(13) 
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becomes 


(  Ou 

po  V  ~dt  +  ^ U ^ =  p°(1~  «voAT)  g  -  Vp. 


(14) 


Applying  the  curl  operator  yields  the  vorticity  equation: 


duj 


—  +  (Vu>)  u  =  (V«)  oj  +  av  o  jVx  (AT  k ) , 


(15) 


where  k  is  the  unit  vector  in  the  vertical  direction,  and  g  =  —  g  k  is  the  gravity 
vector.  The  baroclinic  vorticity  generation  is  now  just  a  function  of  the  horizontal 
temperature  gradient. 

Equations  15  and  10  are  the  coupled  momentum  and  energy  equations  that  are  the 
focus  of  this  work.  Written  in  nondimensional  they  are 


du) 

—  +  (Vo>)  u  =  (Vm)  u>  +  Ri  V  x  (0  k) 

i  +  «.vr=iv>e, 


Of  Pe 

where  0  is  the  nondimensional  temperature 

T  -  Tref 


0  = 


Th.  -  T, 


ref 


Pe  is  the  Peclet  number 


and  Ri  is  the  Richardson  number 


PP  = 


LU 


a 


(16) 


(17) 


(18) 


Ri  g  Qty o  ( Th  ^re/) 


Ur 


(19) 


Th  is  the  high  temperature,  Tri  f  is  the  reference  or  ambient  temperature,  L  is  the 
length  scale,  U  is  the  characteristic  velocity,  and  a  is  the  thermal  diffusivity. 


The  temperature  equation  is  solved  using  a  distribution  of  energy  particles  that  rep¬ 
resent  the  internal  energy  of  the  system.  Since  the  VFM  is  a  Lagrangian  method, 
solving  the  energy  equations  in  this  manner  is  a  natural  choice.  The  velocity  field 
induced  by  the  filaments  convects  the  particles  and  a  probabilistic  approach  is  used 
to  model  the  diffusion  process. 


For  a  given  sensing  volume  Vs,  the  average  temperature  in  that  volume  is  a  function 
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of  the  number  of  energy  particles  in  the  volume  and  the  energy  each  carries.  To  be 
more  precise 

Yf  Ei 

T  =  Tref  +  *,  (20) 

pcvVs 

where  N  is  the  number  of  particles  in  Vs,  Ei  is  the  energy  carried  by  the  ith  particle, 
p  is  the  fluid  density,  cv  is  the  specific  heat  at  constant  volume,  and  Tre/  is  the 
ambient  temperature.  The  relationship  between  the  nondimensional  temperature  0 
and  the  energy  particle  distribution  becomes 


Tr-p  T7i 

=  ~  Eef  =  2^j  Ei 

Th  —  Tref  pcvVsAT ’ 


where  AT  =  T/(  —  Tref. 


(21) 


Assume  that  Et  =  E0  \/Et  e  Vs,  then 


©  = 


NEo 

pcvVsAT' 


(22) 


In  this  case,  the  particle  density  pp  =  { r  that  is  necessary  to  raise  0  by  1 ,  or  T  by 
A Tw,  is 

pcvA  Tw 
Pp =  E  ' 

£/n 


(23) 


Substituting  this  into  the  previous  equations  gives  a  convenient  computational  form 
for  computing  temperature  of  a  distribution  of  particles: 


PpVa  ’ 


0  = 


(24) 


where  E*  =  J|. 


The  resolution  of  the  temperature  and  temperature  gradient  computations  is  a  func¬ 
tion  of  the  size  of  the  sensing  volume  Vs  and  the  number  of  particles  in  Vs  that 
represents  a  nondimensional  temperature  rise  of  1,  that  is  pp  Vs.  These  are  key  pa¬ 
rameters  in  the  numerical  algorithm. 
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3.  Numerical  Algorithm 


3.1  Overview 

Algorithm  1  presents  a  high-level  overview  of  the  vortex  filament  method.  In  step 
1,  the  flow  domain  is  initialized  as  appropriate  with  vorticity  in  the  form  of  fila¬ 
ments  and  temperature  in  the  form  of  energy  particles.  At  each  time  step,  temper¬ 
ature  gradients  are  computed  and  the  associated  vorticity  is  released  into  the  flow 
as  filaments.  Tube  end  points  and  energy  particles  are  then  advected  using  a  fourth- 
order  Runge-Kutta  method  followed  by  a  Monte  Carlo  method  that  diffuses  the 
energy  particles.  Tubes  exceeding  the  length  constraint  are  split,  and  those  that  do 
not  conform  with  the  minimum  length  are  removed  and  the  associated  filament  is 
reconnected.  Finally,  loops  are  detected  and  removed. 

Algorithm  1  Vortex  Filament  Method 
1:  Initialize  vorticity  and  temperature 
2.  for  i  1  y  V steps  tit) 

3:  Compute  temperature  gradients 

4:  Release  filaments 

5:  Advect  tube  end  points  and  energy  particles 

6:  Diffuse  energy  particles 

7:  Split  tubes  that  exceed  the  maximum  length  constraint 

8:  Remove  tubes  that  do  not  conform  to  the  minimum  length  constraint 

9:  Detect  and  remove  loops 


This  algorithm  was  implemented  in  C++  using  a  hybrid  Message  Passing  Interace 
(MPI)  and  OpenMP  parallelization  strategy.  The  remainder  of  this  section  discusses 
the  implementation  details. 

3.2  Data  Distribution 

Filament  operations  and  velocity  calculations  have  competing  requirements  relative 
to  data  distribution  and  parallelization.  Having  all  tubes  in  a  filament  on  the  same 
processor  is  most  efficient  when  searching  for  and  removing  loops,  but  then  large 
filaments  are  likely  spread  over  a  significant  portion  of  the  physical  domain.  On  the 
other  hand,  efficient  velocity  calculations  prefer  that  the  physical  domain  is  decom¬ 
posed  into  small,  contiguous  regions  with  minimal  shared  boundaries.  This  likely 
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leads  to  large  filaments  spread  across  multiple  processors.  There  was  no  attempt  to 
resolve  these  competing  interests  because  a  relatively  simple  algorithm  for  detect¬ 
ing  loops  becomes  complex  if  a  filament’s  tubes  are  widely  distributed;  therefore, 
filaments  are  distributed  over  the  processors  intact.  The  resulting  communication 
cost  from  this  decision  is  minimized  in  some  situations  by  optionally  redistributing 
the  filaments  based  on  the  distribution  of  its  tubes  during  the  previous  velocity  cal¬ 
culation:  the  processor  with  the  most  tubes  gets  the  entire  filament.  The  actual  data 
distribution  for  the  velocity  and  temperature  calculations  depends  on  the  octrees 
and  load  balancing  schemes,  which  are  discussed  next. 

3.3  Octree _ 

Velocity,  temperature,  and  temperature  gradient  calculations  all  use  an  octree  data 
structure  so  a  single  template  framework  was  written  and  used  for  each.  Given 
some  distributed  data  set,  a  global  octree  is  constructed  by  distributing  the  data 
across  multiple  MPI  processes  and  building  local  octrees  that  share  global  proper¬ 
ties.  The  data  are  distributed  by  choosing  some  level  of  the  global  octree  to  be  the 
load-balance  level  and,  based  on  workload,  assigning  the  tree  nodes  at  that  level, 
and  consequently  the  subtrees  rooted  at  those  nodes,  to  the  processors.  Two  load¬ 
balancing  algorithms  are  available  for  this:  one  based  on  distributing  the  subtrees 
in  Morton  Key  order  to  maintain  contiguous  data  sets  on  each  processor,  and  an¬ 
other  that  simply  sorts  the  nodes  by  the  amount  of  work  and  distributes  the  nodes  in 
round-robin  fashion.  Once  built,  octree  controller  classes  manage  the  local  octrees 
and  the  communication  between  them  as  required  for  the  specific  calculation.  Tem¬ 
plate  arguments  define  the  octree  data  types  allowing  the  framework  to  be  generic. 

Construction  of  the  global  octree  is  similar  for  each  specific  calculation.  Aside  from 
differing  data  types,  some  minor  differences  might  include  limits  to  the  maximum 
tree  depth  and  the  load  balancing  algorithm  used.  Before  building  the  octree,  all 
processors  must  agree  on  global  data  bounds,  the  load-balance  level,  and  a  common 
tree  structure  above  the  load-balance  level.  (The  last  requirement  is  only  necessary 
when  used  in  the  fast  multipole  method.  It  makes  it  convenient  for  exchanging  mul¬ 
tipole  expansions.)  Data  on  each  processor  are  sorted  according  to  its  Morton  key 
based  on  the  global  data  bounds  and  the  lowest  possible  level  of  the  tree.  Each  local 
data  point  is  assigned  to  a  load-balance  node  and  metadata  describing  the  local  data 
distribution  are  exchanged  with  all  processors.  The  chosen  load  balancing  algorithm 
then  uses  these  metadata  to  decide  on  an  optimal  distribution  and  redistributes  the 
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data  accordingly.  Once  redistributed,  each  processor  constructs  a  local  octree.  This 
is  summarized  in  Algorithm  2. 


Algorithm  2  Global  Octree  Construction 

Assumption:  A  nonoptimal  data  distribution  exists 
1:  Compute  local  data  bounds 
2:  Exchange  local  bounds  and  set  global  data  bounds 
3:  Set  the  load-balance  level 

4:  Locally  sort  and  assign  data  to  a  load-balance  node 
5:  Exchange  load  balance  meta  data 
6:  Invoke  load  balance  algorithm 

7:  Redistribute  the  data  based  on  output  from  the  load  balance  algorithm 
8:  Build  local  octrees 


This  octree  framework  is  used  in  the  FMM  when  computing  velocities  during  the 
advection  step  and  when  computing  temperature  and  temperature  gradients  as  part 
of  the  filament  release  process. 

3.4  Advection _ 

At  each  time  step,  tube  end  points  and  energy  particles  are  advected  using  a  fourth- 
order  Runge-Kutta  routine  that  requires  4  velocity  calculations  per  time  step.  If  M 
is  the  number  of  vortex  tubes  and  N  is  the  number  of  unique  tube  end  points  plus 
energy  particles,  then  the  velocity  calculation  is  0{MN).  This  is  cost  prohibitive 
for  any  reasonable  size  problem  so  the  FMM6  is  used  to  reduce  the  velocity  com¬ 
putation  to  0(A). 

The  intent  of  this  work  was  not  to  develop  a  new  FMM  since  much  work  had  al¬ 
ready  been  done  in  this  area.  The  ExaFMM,1  circa  September  2016,  was  chosen 
for  this  and  integrated  into  the  VFM.  Some  difficulties,  many  associated  with  non- 
uniformly  distributed  data  as  found  in  most  problems  of  interest,  were  encountered 
so  a  new  FMM  was  written  based  in  part  on  the  ExaFMM.  This  new  FMM  uses 
the  ExaFMM  spherical  harmonics  routines  and  mimics  its  traversal  algorithms.  It 
is  written  as  a  template  C++  class  and  uses  the  octree  framework  discussed  pre¬ 
viously.  Template  arguments  define  the  kernel  type,  for  example,  Biot-Savart  and 
Laplace,  and  the  associated  data  types.  It  is  a  list-based  dual  tree  approach  where 
separate  octrees  are  constructed  for  the  sources  and  the  targets.  The  dual  tree  ap¬ 
proach  is  well  suited  for  the  VFM  because  the  depth  of  the  source  tree  is  limited  by 
the  smoothing  function  but  not  so  for  the  target  tree.  When  N  '^>  M,  the  hierarchy 
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of  the  target  tree  can  be  exploited  when  translating  the  local  expansions  to  the  bot¬ 
tom  of  the  target  octree,  which  in  general  is  at  a  much  lower  depth.  The  details  of 
this  FMM  will  be  documented  in  a  separate  technical  report. 

3.5  Temperature  Gradient 

The  temperature  and  temperature  gradient  calculations  start  by  constructing  an  oc¬ 
tree  that  captures  all  the  energy  particles  and  whose  depth  is  limited  by  the  size 
of  the  input  sensing  volume.  Figure  4a  shows  an  example  of  a  spherical  cloud  of 
energy  particles  enclosed  by  the  sensing  volume  octree.  Global  bounds  of  the  oc¬ 
tree  are  chosen  so  that  node  volumes  at  the  lowest  level  of  the  tree  are  equal  to  the 
desired  sensing  volume.  An  input  parameter,  ntinf,  defines  a  minimum  number  of 
energy  particles  above  which  the  temperature  is  considered  greater  than  zero.  This 
controls  the  growth  of  the  tree  such  that  the  number  of  energy  particles  in  nodes  at 
the  sensing  volume  level  will  always  contain  greater  than  or  equal  to  this  amount 
or  will  have  at  least  one  sibling  that  does.  Once  constructed,  it  is  a  simple  matter  to 
sum  the  energy  of  all  the  particles  in  a  sensing  volume  to  compute  the  temperature. 


(a)  Sensing  volume  octree 


(b)  Least  squares  domain 


(c)  Filament  release 


Fig.  4  Temperature  and  baroclinic  filament  release  calculations  utilize  an  octree  with  leaf  nodes 
at  the  sensing  volume  (4a),  a  least  squares  curve  fit  for  temperature  using  the  node  (red)  and 
its  nearest  neighbors  (4b),  and  a  filament  released  inside  a  sensing  volume  (or  release  volume) 
(4c). 


Temperature  gradients  are  computed  at  the  center  of  each  sensing  volume  using  the 
gradient  of  a  full  second-order  polynomial  approximation  of  temperature.  Figure  4b 
shows  an  example  where  the  red  node,  surrounded  by  its  nearest  neighbors,  is  where 
the  temperature  gradient  is  to  be  computed.  The  coefficients  of  the  polynomial  are 
generated  using  a  least  squares  curve  fit  with  data  from  the  node  and  its  nearest 
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neighbors.  If  one  of  its  26  near  neighbors  does  not  exist  in  the  octree  it  is  because 
the  nondimensional  temperature  there  is  0. 


3.6  Filament  Release _ 

Once  temperature  gradients  are  computed  the  associated  vorticity  is  know  and  can 
be  released  into  the  flow  as  filaments.  This  can  be  done  by  constructing  filaments 
inside  the  sensing  volume  itself,  or  by  partitioning  the  domain  into  equivalent  vol¬ 
umes  in  some  other  way  and  releasing  from  there.  It  was  found  that  when  a  problem, 
such  as  a  spherical  thermal  bubble,  has  some  symmetry  not  aligned  with  axes,  using 
sensing  volumes  that  are  aligned  with  the  axes  to  construct  the  filaments  results  in 
solutions  that  show  some  grid  dependency.  This  can  be  mitigated  by  either  reducing 
the  sensing  volume  further,  which  can  significantly  increase  the  number  of  energy 
particles  needed,  or  by  constructing  release  volumes  that  reflect  the  symmetry  of  the 
problem.  Temperature  gradients  can  be  computed  at  the  center  of  any  chosen  release 
volume  using  the  appropriate  polynomial  approximation  for  temperature.  Release 
volumes  that  have  radial  symmetry  were  used  for  the  spherical  bubble  calculations 
presented  in  the  Section  4. 

Filaments  are  constructed  by  intersecting  the  release  volumes  with  a  line  that  passes 
through  the  center  and  is  parallel  to  the  vorticity  vector:  see,  for  example,  Fig.  4c. 
The  magnitude  of  the  vorticity,  the  length  of  intersecting  line  and  the  release  volume 
determine  the  circulation  of  the  new  filament.  If  the  new  filament  is  longer  than  the 
maximum  tube  length  it  is  subdivided.  For  release  volumes  that  are  not  cubes,  the 
same  procedure  is  applied  to  an  equivalent  cube  centered  over  the  release  volume. 

3.7  Loop  Removal 

Loop  removal  models  small-scale  dissipation  while  limiting  the  growth  of  the  num¬ 
ber  of  tubes;  thus,  it  is  an  essential  part  of  the  method.  The  tube  end  points  of  a 
filament  are  used  for  loop  detection.  When  the  distance  between  any  2  of  these 
points  is  less  than  the  distance  criterion  di,  that  portion  of  the  filament  is  considered 
a  loop,  the  tubes  between  them  are  removed,  and  the  filament  is  reconnected.  The 
distance  criterion  di  is  based  on  the  total  length  l  of  the  potential  loop  and  an  input 
parameter  a,  and  is  computed  by 

di  = - l,  0  <  a  <  1.  (25) 

1  —  a 
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For  small  a,  di  is  approximately  al.  A  typical  value  is  0.0025.  An  additional  pa¬ 
rameter,  ls,  controls  detection  of  small  loops.  If  l  <  ls  then  di  =  1/2.  Typically,  ls  is 
very  small  if  used  at  all.  When  all  loops  are  detected  and  removed,  if  the  number  of 
tubes  remaining  is  less  than  a  given  percentage  of  the  original  number  of  tubes  the 
entire  filament  is  removed.  This  is  controlled  by  input  parameter  /3. 

The  loop  removal  algorithm  as  described  is  shown  in  Algorithm  3.  The  input  pa¬ 
rameters  are  n0,  na,  ls,  a,  and  f3 ,  where  nQ  is  the  minimum  number  of  tubes  the 
filament  must  have  before  being  checked  for  loops  and  na  is  the  minimum  number 
of  tubes  that  any  loop  must  possess.  The  other  input  parameters  are  as  described 
previously. 


Algorithm  3  Loop  Removal 


Input:  na,ls,a,/3 


1 

2: 

3: 

4: 

5: 

6 

7: 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22: 
23: 
24 
25: 
26 


for  each  filament  do 

if  Number  of  tubes  >  n0  then 

ne  =  number  of  filament  end  points;  i  =  0 
while  i  <  ne  —  na  do 

j  =  i  +  na;  l  =  YjJn=i  found  =  false 
while  not  found  and  j  <  ne  do 

1  =  1  +  lj-i‘,  d  =  |xj  —  Xj| 

if  l  <  L  then 


d  —  1 

U.J  —  2 


else 

di  =  t —l 

1  1 —a 

if  d  <  di  then 

flag  tube  end  points  between  i  and  j  for  removal 
ifi==0ori  —  1  was  marked  for  removal  then 
mark  i  for  removal 

if  j  ==  ne  —  lori  was  marked  for  removal  then 
mark  j  for  removal 
i  =  j,  found  =  true 
else 

j  =  3  +  1 

if  found  =  false  then 
i  =  i  +  1 

if  fraction  of  end  points  marked  for  removal  >  (3  then 
remove  the  entire  filament 
else 

remove  end  points  marked  for  removal  and  reconnect 
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4.  Results 


VFM  simulations  of  the  7  ellipsoidal  bubbles  reported  in  Shapiro7  were  performed. 
The  initial  temperature  distribution  in  each  bubble  was  defined  by  Eqs.  3.4  through 
3.6  of  Shapiro7  and  are  repeated  here: 


T 


B 


[3 

LI 


-SObMS-1 


exp 


(26) 


where  Lx,  Ly  and  Lz  are  the  semi-principal  axes  of  the  ellipsoid, 


x 2  y2 

M  +  M 


(27) 


and 


B  = 


A  Tb 

l  +  Ll/Ll  +  Ll/Ll' 


(28) 


Table  1  defines  the  semi-principal  axes  and  A Tb;  and  the  maximum  temperature 
rise  for  each  case. 


Each  problem  was  nondimensionalized  using  a  reference  length  Lref  defined  by 
Lref  =  max  ( Lx,Ly,Lz )  and  a  reference  time  tref  given  by  the  time  it  would  take 
for  a  bubble  at  constant  temperature  to  rise  a  distance  Lref  if  the  acceleration  due 
to  buoyancy  remained  constant.  The  acceleration8  used  was 


_  mb-md  _  A Tb 

9  .  9  r>  rri  .  A  rji  5 

mb  +  md  2  +  A  Tb 


(29) 


where  mb  is  the  mass  of  the  bubble  and  rnd  is  mass  of  the  displaced  fluid.  The  ref¬ 
erence  time  is  then  tref  =  \JC2  Lr(,f /ab,  the  reference  velocity  is  Uref  =  Lref/tref, 
and  the  reference  temperature,  Tref,  is  30076  for  all  cases.  The  Richardson  number 
written  in  terms  of  Tre/  and  A Tb  is 


Ri 


2-^ 

-L  ref 


4  + 


2ATb 

Tref 


Table  1  lists  the  reference  quantities  and  Richardson  number  for  each  case. 
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Table  1  Parameters  and  reference  quantities  for  the  Shapiro-Kanak  bubbles.  Tref  =  300  K 
for  all  cases 


Case 

Lx  (m) 

Ly  (m) 

Lz  (m) 

A Tb  (K) 

Lref  (m) 

tref  00 

Uref  (m/s) 

Ri 

CNTRL 

64.24 

64.24 

64.24 

1.5 

64.24 

72.5 

0.886 

4.01 

EXPT1 

26.84 

64.24 

64.24 

1.5 

64.24 

72.5 

0.886 

4.01 

EXPT2 

20.84 

49.88 

64.24 

1.5 

64.24 

72.5 

0.886 

4.01 

EXPT3 

20.84 

64.24 

49.88 

1.5 

64.24 

72.5 

0.886 

4.01 

EXPT4 

26.84 

64.24 

26.84 

1.5 

64.24 

72.5 

0.886 

4.01 

EXPT5 

20.84 

49.88 

19.23 

1.5 

49.88 

63.89 

0.781 

4.01 

EXPT6 

20.84 

49.88 

19.23 

3.0 

49.88 

45.23 

1.103 

4.02 

VFM  algorithm  parameters  used  are  given  in  Tables  2  through  4.  Uniform  energy 
particles  were  used  for  all  cases;  that  is,  the  energy  contained  in  each  particle  is 
the  same.  The  parameters  in  Table  3  imply  that  8,192  particles  are  required  to  raise 
the  nondimensional  temperature  in  a  sensing  volume  by  1 .  For  all  cases  considered, 
convection  dominated  and  diffusion  was  turned  off. 

The  thermal  bubbles  were  initialized  with  a  distribution  of  energy  particles  that 
represented  the  temperature  given  by  Eq.  26.  This  was  done  by  covering  the  bubble 
with  a  net  of  cubes,  each  of  whose  volume  was  ^  °f  the  sensing  volume.  The 
temperature  in  each  cube  was  considered  constant  and  the  appropriated  number  of 
energy  particles  was  deposited  randomly  into  the  cube.  Figure  5  shows  3  views  of 
the  initial  cloud  of  particles  and  Table  5  gives  the  total  number  of  energy  particles 
needed  for  each  case. 

Each  case  was  run  to  an  actual  time  of  216  s  using  a  constant  nondimensional  time 
step  dt  =  0.033103.  The  reference  times  for  cases  EXPT5  and  EXPT6  implied  that 
102  and  144  iterations,  respectively,  were  needed  to  reach  the  simulation  time.  The 
other  cases  required  90  iterations.  As  can  be  seen  from  Table  1,  doubling  A Tb  had 
a  minimal  effect  on  the  Richardson  number  and  a  large  effect  on  the  reference  time. 
Because  of  this,  EXPT5  and  EXPT6  are  essentially  the  same  case  run  to  different 
nondimensional  times.  Figure  6  shows  how  the  initial  cloud  of  energy  particles 
evolved  after  216  s. 
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Table  2  Vortex  tube  parameters 


Parameter 

Value 

Pc 

0.01 

Imax 

0.025 

Imin 

0.0001 

Table  3  Energy  particle  parameters 


Parameter 

Value 

Pp 

16  x  10° 

Vs 

512  x  10-6 

Tltinf 

8 

Table  4  Loop  removal  parameters 


Parameter 

Value 

n0 

5 

na 

5 

h 

0.0 

a 

0.005 

P 

0.8 

Table  5  Number  of  energy  particles  used 


Case 

Number 

CNTRL 

89,322,424 

EXPT1 

48,777,768 

EXPT2 

31,832,360 

EXPT3 

31,832,360 

EXPT4 

18,383,048 

EXPT5 

17,212,408 

EXPT6 

17,212,408 
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EXPT1 


EXPT2 


EXPT3 


EXPT4 


EXPT5 


EXPT6 


(a) 


I 


(b) 


i 

I 

I 

I 

I 


Fig.  5  Initial  cloud  of  energy  particles  as  viewed  in  the  (a)  —  x,  (b)  +y,  and  (c)  — z  directions 
for  the  1  spherical  and  6  ellipsoidal  thermal  bubbles 
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EXPT3 


EXPT4 


EXPT5 


EXPT6  JL 


(a) 


(b) 


Fig.  6  Same  as  Fig.  5  but  at  t  =  216  s 
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Figures  7  through  13  show  the  distribution  of  vortex  filaments  and  tubes  for  each 
case  at  the  final  time.  In  Figs.  8  through  13,  subfigure  (a)  shows  the  filaments  as 
viewed  in  the  —  x  direction  and  subfigure  (b)  shows  them  as  viewed  in  the  +y  di¬ 
rection.  The  filaments  are  color  coded  with  circulation  using  the  log  color  scale 
shown  in  Fig.  7b.  The  strong  filaments  are  predominately  at  the  top  of  each  bubble 
with  the  strongest  ones  in  the  interior,  which  cannot  be  seen  in  these  figures.  Sub¬ 
figures  (c)  and  (d)  of  Figs.  7  through  13  show  how  circulation  is  distributed  with 
respect  to  the  filaments  and  tubes,  respectively,  and  (e)  shows  the  distribution  of  the 
number  of  tubes  per  filament. 
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(a) 


circulation 

B-1e-02 


1e-03 


1e-04 


(b) 


(c)  (d)  (e) 

Fig.  7  CNTRL  filaments  at  t  =  216  s.  (a)  113,085,992  tubes  colored  by  circulation  and  viewed 
in  the  — x  direction,  (b)  Log  scale  color  map  for  the  circulation,  (c),  (d),  and  (e)  Distributions 
of  circulation  and  tubes  per  filament. 


(a) 


(b) 


(c)  (d)  (e) 

Fig.  8  EXPT1  filaments  at  t  =  216  s.  (a)  334,298,530  tubes  colored  by  circulation  and  viewed  in 
the  —x  direction,  (b)  Filaments  as  viewed  in  the  +y  direction,  (c),  (d),  and  (e)  Distributions  of 
circulation  and  tubes  per  filament. 
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(a) 


(b) 


(c) 


(d) 


(e) 


Fig.  9  Similar  to  Fig.  8  for  EXPT2  with  364,379,492  tubes 


(a) 


(b) 


10°  101  102  103  104 
tubes  per  filament 


(C) 


(d) 


(e) 


Fig.  10  Similar  to  Fig.  8  for  EXPT3  with  323,079,119  tubes 
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(a) 


tubes  per  filament 


(C) 


(d) 


(e) 


Fig.  11  Similar  to  Fig.  8  for  EXPT4  with  77,288,312  tubes 


(a) 


(b) 


(c) 


(d) 


(e) 


Fig.  12  Similar  to  Fig.  8  for  EXPT5  with  107,678,960  tubes 
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(a) 


(b) 


10°  101  102  103  104 
tubes  per  filament 


(C) 


(d) 


(e) 


Fig.  13  Similar  to  Fig.  8  for  EXPT6  with  533,352,010  tubes 


The  lower  limit  in  the  circulation  distribution  charts  reflect  the  limit  imposed  on 
weak  filaments.  Any  filament  generated  with  circulation  less  than  6.0  x  10-5  was 
discarded.  Ad  hoc  calculations  indicated  that  this  cutoff  level  resulted  in  very  few 
discarded  filaments  with  little  effect  on  the  solution;  however,  this  needs  to  be  ex¬ 
amined  further.  The  drop  off  at  the  high  end  is  a  measure  of  the  extent  of  the  largest 
temperature  gradients  encountered.  The  maximum  circulation  for  any  case  is  a  func¬ 
tion  of  the  temperature  gradient,  release  volume,  and  time  step.  The  release  volumes 
were  on  the  order  of  the  sensing  volume  and  were  similar  across  all  cases,  as  was 
the  time  step.  The  drop  off  in  the  distribution  of  the  number  of  tubes  per  filament, 
subfigure  (e),  is  due  to  loop  removal.  Larger  filaments  have  a  higher  probability  of 
forming  loops. 

Figure  14  shows  how  the  number  of  filaments,  tubes  and  tubes  removed  from  loops 
increased  during  the  simulations.  At  the  final  time  step,  tube  growth  remained  ex¬ 
ponential.  More  work  needs  to  be  done  to  limit  this  growth  through  aggressive  loop 
removal,  merging  of  tubes,  or  other  methods,  so  that  these  calculations  can  be  ad¬ 
vanced  to  a  fully  turbulent  state  with  more  reasonable  tube  growth. 
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(a)  filaments 


(b)  tubes 


(c)  loop  removal 


Fig.  14  VFM  calculation  statistics  vs.  iteration  number:  (a)  is  the  number  of  filaments,  (b)  is 
the  number  of  vortex  tubes,  and  (c)  is  the  number  of  tubes  removed  from  loops 
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To  compare  the  maximum  vertical  velocity  and  vorticity  predicted  by  the  VFM  to 
the  same  presented  in  Figs.  5  and  7  in  Shapiro,7  a  uniform  3-D  grid  of  size  121  x 
121  x  201  with  an  nondimensional  grid  size  of  0.05  was  constructed  around  each 
bubble.  The  velocity  and  vorticity  were  then  computed  at  each  grid  point  and  the 
maximum  of  the  vertical  component  of  each  was  found.  The  comparisons  are  shown 
in  Fig.  15  in  dimensional  quantities.  The  black  vertical  axis  on  the  left  is  the  velocity 
scale  and  the  red  vertical  axis  on  the  right  is  the  vorticity  scale.  The  curves  are 
similarly  color  coded.  The  solid  lines  are  VFM  simulations  and  the  dashed  lines  are 
numbers  picked  off  Figs.  5  and  7  in  Shapiro.7  The  velocity  compares  well  with  the 
published  results  with  excellent  agreement  on  the  spherical  bubble.  At  early  times, 
the  maximum  vertical  vorticity  is  always  larger  than  Shapiro.7  This  is  most  likely 
due  to  the  nondiffusive  nature  of  the  VFM  method  that  allows  tubes  to  reorient  even 
with  minimal  perturbation  as  is  the  case  for  early  times.  At  later  times,  these  curves 
tend  to  come  closer  together. 

Figure  16  shows  the  height  reached  by  the  bubbles  as  a  function  of  time.  The  curves 
with  "_sk"  appended  to  their  labels  were  generated  by  estimating  the  height  from 
Fig.  2  of  Shapiro7  for  the  CNTRL  case  and  Fig.  3  for  the  EXPT1  case.  Excellent 
agreement  with  these  2  cases  is  shown. 

Figure  17  shows  the  nondimensional  temperature  difference  in  the  at  x  =  0  plane 
for  the  CNTRL  bubble  at  nonscaled  times  t  =  96  through  t  =  216  in  increments  of 
24  s.  Figures  18  through  23  show  the  same  for  the  nonspherical  bubbles  in  addition 
to  the  temperature  difference  in  the  y  =  0  plane.  Figures  17  and  18  can  be  compared 
with  Figs.  2  and  3  in  Shaprio.7 
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(a)  CNTRL 


(b)  EXPT1 


(c)  EXPT2 


(d)  EXPT3 


(e)  EXPT4 


(f)  EXPT5 


(g)  EXPT6 


Fig.  15  Maximum  vertical  velocity  and  vorticity. 
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Fig.  16  Height  of  the  thermal  bubbles  vs.  time 
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(a)  t  =  96 


(b)  t  =  120  (c)  t  =  144 


(d)  t  =  168  (e)  t  =  192  (f>  t  =  216 


Fig.  17  CNTRL  bubble  temperature  rise  in  the  x  =  0  plane.  See  Fig.  24  for  the  color  scale. 
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(g)  t  =  96  (h)  t  =  120  (i)  t  =  144  (j)  t  =  168  (k)  t  =  192 


(1)  t  =  216 


Fig.  18  EXPT1  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 
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(g)  t  =  96  (h)  t  =  120  (i)  t  =  144  (j)  t  =  168  (k)  t  =  192  (1)  t  =  216 


Fig.  19  EXPT2  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 
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(g)  t  =  96  (h)  t  =  120  (i)  t  =  144  (j)  t  =  168  (k)  t  =  192  (1)  t  =  216 


Fig.  20  EXPT3  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 
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(g)  t  =  96  (h)  t  =  120  (i)  t  =  144  (j)  t  =  168  (k)  t  =  192  (1)  t  =  216 


Fig.  21  EXPT4  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 
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Fig.  22  EXPT5  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 


Approved  for  public  release;  distribution  is  unlimited. 


36 


0  50  100 

y  (m) 


(a)  t  =  96 


(b)  t  =  120  (c)  t  =  144  (d)  t  =  168  (e)  t  =  192 


0  50  100 

y  (m) 

(f)  t  =  216 


0  50  100 


0  50  100 


0  50  100 


0  50  100 


0  50  100 


0  50  100 


(§)  t  =  96 


(h)  t  =  120 


(i)  t  =  144 


(j) t  =  168  (k) t  =  192 


(1)  t  =  216 


Fig.  23  EXPT6  bubble  temperature  rise  in  the  x  =  0  plane,  (a)-(f),  and  the  y  =  0  plane,  (g)-(l). 
See  Fig.  24  for  the  color  scale. 
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Fig.  24  Nondimensional  temperature  scale  for  Figs.  17-23 
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5.  Conclusions  and  Path  Forward 


The  VFM  offers  a  methodology  for  solving  the  heretofore  intractable  problem  of 
capturing  the  physics  of  nondiffusive  energy  transfer  and  vortex  dynamics  in  high 
Re  turbulent  flows.  Couple  with  the  energy  equation,  the  results  presented  here 
show  that  the  method  is  also  able  to  capture  the  correct  physics  of  thermally  driven 
flows.  The  vertical  velocity,  structural  form,  and  temperature  distribution  of  the 
thermal  bubbles  match  published  calculations.  This  new  method  has  the  potential 
for  accurately  predicting  turbulent  atmospheric  flow  environments  in  ways  not  seen 
before  and  provide  the  means  for  better  understanding  such  flows. 

A  separate  validation  study  of  the  method,  not  reported  here,  focused  on  simulations 
of  isotropic  turbulence  in  a  periodic  box.  Results  show  that  the  method  produces 
isotropic  turbulence,  for  example,  2-point  velocity  correlations  obey  isotropic  laws. 
These  results  will  be  published  in  a  separate  report. 

The  capability  to  simulate  heated  surfaces  has  been  added,  and  future  work  will 
include  simulations  of  thermally  driven  flows  initialized  from  vertical  and  horizon¬ 
tal  heated  surfaces.  The  latter  will  provide  results  that  can  be  compared  with  data 
collected  at  Meteorological  Sensor  Array  at  the  White  Sands  Missile  Range  in  New 
Mexico.  Additional  planned  capabilities  include  adding  viscous  surfaces  for  simu¬ 
lations  with  complex  geometries,  2-way  coupled  particles  for  simulations  with  dust 
and  chemical  agents,  and  the  addition  of  a  3-D  radiative  transfer  scheme. 

Further  research  is  needed  to  contain  the  growth  of  vortex  tubes  through  refinement 
of  the  loop  removal  process  or  other  means  such  as  combining  tubes.  The  effects 
on  the  underlying  physics  when  energy  is  lost  due  to  loop  removal  needs  further 
study. 
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Appendix.  Integrals 
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The  exact  integration  of  the  first  2  integrals  for  uz  in  Eq.  9  are 


(A-l) 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3-D 

3-dimensional 

ABL 

atmospheric  boundary  layer 

DNS 

direct  numerical  simulation 

FMM 

Fast  Multipole  Method 

GPUs 

graphic  processing  units 

LES 

Large  Eddy  Simulation 

M-0 

Monin-  Obukhov 

MPI 

Message  Passing  Interface 

Re 

Reynolds  number 

SFS 

subfilter  scale 

SGS 

subgrid  scale 

VFM 

Vortex  Filament  Method 
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